home *** CD-ROM | disk | FTP | other *** search
/ Best of Shareware / Best of PC Windows Shareware 1.0 - Wayzata Technology (7111) (1993).iso / mac / DOS / PROGRAMG / GRAD / FUNC4.FOR < prev    next >
Text File  |  1988-08-09  |  3KB  |  99 lines
  1.       SUBROUTINE VALUE(A,B,C,A44,C44,FN,INF)
  2.       IMPLICIT REAL*8(A-H,O-Z)
  3.       REAL*4 X
  4.       DIMENSION A(2,4),B(2),C(2,4)
  5. C ****** SUBST. MAX. NUMBER OF OBSERVATIONS ******
  6.       COMMON X(13,1800),NOBS,B3
  7.       COMMON/TRAPCODE/ITR
  8.       COMMON /VRNCS/ W1,W2,W3
  9.       ITR=0
  10. C -- COMPUTE EIGENVALUES AND EIGENVECTORS -
  11.       DET=A(1,1)*A(2,2)-A(1,2)*A(2,1)
  12.       HFTRA=0.5*(A(1,1)+A(2,2))
  13.       DISCR=HFTRA**2-DET
  14.       IF(DISCR.GT.0.) GO TO 1
  15.          WRITE(3,11)
  16.    11    FORMAT(' COMPLEX EIGENVALUES')
  17.          INF=1
  18.          RETURN
  19.     1 IF(HFTRA.lt.0.) goto 91
  20.       XL1=HFTRA+DSQRT(DISCR)
  21.  91   continue      
  22.       IF(HFTRA.ge.0.) goto 92
  23.       XL1=HFTRA-DSQRT(DISCR)
  24.  92   continue
  25.       XL2=DET/XL1
  26.       P11=XL1-A(2,2)
  27.       P12=A(1,2)
  28.       P13=(A(1,3)*P11+A(2,3)*P12)/(XL1+1.)
  29.       P14=(A(1,4)*P11+A(2,4)*P12)/(XL1-A44)
  30.       P21=A(2,1)
  31.       P22=XL2-A(1,1)
  32.       P23=(A(1,3)*P21+A(2,3)*P22)/(XL2+1.)
  33.       P24=(A(1,4)*P21+A(2,4)*P22)/(XL2-A44)
  34. C -- COMPUTE LOG-LIKELIHOOD -
  35.       DETC=C(1,1)*C(2,2)-C(1,2)*C(2,1)
  36.       IF(ABS(DETC).GT.0.) GO TO 3
  37.          WRITE(3,33)
  38.    33    FORMAT(' C IS SINGULAR')
  39.          INF=1
  40.          RETURN
  41.     3 DETP=P11*P22-P12*P21
  42.       IF(ABS(DETP).GT.0.) GO TO 4
  43.          WRITE(3,44)
  44.    44    FORMAT(' P IS SINGULAR')
  45.          INF=1
  46.          RETURN
  47.     4 continue
  48.       FN=-DLOG(ABS(DETC*DETP))
  49.       W1=0.
  50.       W2=0.
  51.       W3=0.
  52.       PB1=P11*B(1)+P12*B(2)+P13*B3-P14*A44
  53.       PB2=P21*B(1)+P22*B(2)+P23*B3-P24*A44
  54. C -- LOOP BEGINS -
  55.       DO 10 I=1,NOBS
  56.       IF(ITR.EQ.1) GO TO 2
  57.       Y04=1.0-(1.0-DEXP(X(11,I)*C44))*DEXP(X(12,I)*(-A44))
  58.       Y4=1.0-(1.0-DEXP(X(11,I)*C44))*DEXP(X(13,I)*(-A44))
  59.       ALX04=DLOG(Y04)/C44
  60.       ALX4=DLOG(Y4)/C44
  61.       ALXC1=C(1,1)*X(1,I)+C(1,2)*X(2,I)+C(1,3)*X(3,I)+C(1,4)*ALX4
  62.       ALXC2=C(2,1)*X(1,I)+C(2,2)*X(2,I)+C(2,3)*X(3,I)+C(2,4)*ALX4
  63.       Y1=DEXP(ALXC1)
  64.       Y2=DEXP(ALXC2)
  65.       Y01=DEXP(C(1,1)*X(4,I)+C(1,2)*X(5,I)+C(1,3)*X(6,I)+C(1,4)*ALX04)
  66.       Y02=DEXP(C(2,1)*X(4,I)+C(2,2)*X(5,I)+C(2,3)*X(6,I)+C(2,4)*ALX04)
  67.       EL1T=DEXP(X(7,I)*(-XL1))
  68.       EL2T=DEXP(X(7,I)*(-XL2))
  69.       E1=P11*Y1+P12*Y2+P13*X(8,I)+P14*Y4
  70.      1     -EL1T*(P11*Y01+P12*Y02+P13*X(9,I)+P14*Y04)
  71.       E2=P21*Y1+P22*Y2+P23*X(8,I)+P24*Y4
  72.      1     -EL2T*(P21*Y01+P22*Y02+P23*X(9,I)+P24*Y04)
  73.       IF(ABS(XL1).GT.1.D-22) GO TO 222
  74.          R1=-X(7,I)
  75.          E1=E1-R1*PB1
  76.          GO TO 333
  77.  222  R1=(EL1T*EL1T-1.)/(XL1+XL1)
  78.       E1=E1+(1.-EL1T)/XL1*PB1
  79.  333  IF(ABS(XL2).GT.1.D-22) GO TO 444
  80.          R2=-X(7,I)
  81.          E2=E2-R2*PB2
  82.          GO TO 555
  83.  444  R2=(EL2T*EL2T-1.)/(XL2+XL2)
  84.       E2=E2+(1.-EL2T)/XL2*PB2
  85.  555  R3=(1.-DEXP(2.D0*X(7,I)))*0.5
  86.       W1=W1+E1*E1/R1
  87.       W2=W2+E2*E2/R2
  88.       W3=W3+X(10,I)*X(10,I)/R3
  89.       FN=FN+(0.5*DLOG(R1*R2*R3)-ALXC1-ALXC2)/NOBS
  90.    10 continue
  91. C -- LOOP ENDS -
  92.       FN=FN+0.5*(DLOG(W1*W2*W3))
  93.     2 CONTINUE
  94.       INF=ITR
  95.  
  96.       RETURN
  97.       END
  98.  
  99.